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ABSTRACT 

Aims. We investigate the statistics for the spatial and temporal distribution of the energy input into the corona in a three-dimensional 
magneto-hydrodynamical (3D MHD) model. The model describes the temporal evolution of the corona above an observed active 
region. The model is driven by photospheric granular motions that braid the magnetic field lines. This induces cuiTents that are 
dissipated, thereby leading to transient heating of the coronal plasma. We evaluate the transient heating as subsequent heating events 
and analyze their statistics. The results are then interpreted in the context of observed flare statistics and coronal heating mechanisms. 
Observed solar flares and other smaller transients cover a wide range of energies. The frequency distribution of energies follow a 
power law, the lower end of the distribution given by the detection limit of current instrumentation. One particular heating mechanism 
is based on the occuiTence of so-called nanoflares, i.e. very low-energy deposition events. 

Methods. To conduct the numerical experiment we use a high-order finite-difference code that solves the partial differential equations 
for the conservation of mass, the momentum and energy balance, and the induction equation. The energy balance includes Spitzer 
heat conduction and optically thin radiative losses in the corona. 

Results. The temporal and spatial distribution of the Ohmic heating in the 3D MHD model follows a power law and can therefore be 
understood as a system in a self-organized critical state. The slopes of the power law are similar to the results based on observations of 
flares and smaller transients. We find that the coronal heating is dominated by events similar to the so-called nanoflares with energies 
on the order of 10" J or 10^"* erg. 
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1. Introduction 

In the upper atmosphere of the Sun and stars, energy is released 
in a transient fashion. Strong flares are just the tip of the iceberg, 
where the plasma can be heated up to 10 MK or more and par- 
ticles are accelerated to relativistic speeds. More numerous are 
smaller energy releases that cause smaller brightenings in X-rays 
and extreme UV, down to the current detection limit. There is a 
continuous connection from these smallest observable brighten- 
ings, often called microflares, to the largest of flares in the sense 
that the distribution in energy of all these events roughly fol- 
lows a power law. Theoretical arguments have been proposed 
that even s maller transie nt releases of energy have to exist, the 
nanoflares (|Parkerlll988b . which could play a major role in the 
heating of the corona. In this study we investigate the energy 
distribution of these low energy releases as they are found in 3D 
magnetohydrodynamics (MHD) models. 

Rosner & Vaianal (Il978h have already investigated the tran- 



sient events of the Sun, flaring stars, and other X-ray sources. 
They found that the energy distribution of these events 
roughly follows a pow er law. To explain the flare statistics 
iRosner & Vaianal (1 19781) provided a theoretical model that led to 
a power-law distribution of the flare parameters. This model was 
based on the assumptions that flaring is a stochastic relaxation, 
the waiting time between two events is an independent parame- 
ter, and the energy release during a flare is large in comparison 
to the energy content after the flare. The power-law nature is not 
only observable for X-ray sources, but also solar radio burs ts 
show such a distribution as already studied bv lAkabanel(ll956h . 

These results were followed by numerous observations o f 
solar X-ray flares (ICrosbv et aI.lll993[[T99ilChriste et al.lEoOSi) . 
They have found that not only the flare energy but also other flare 



parameters, such as the flare peak flux or flare duration, follow 
a power law. A review of the statistics o f flares and their sm all 
counterparts, the microflares, is given bv lHannahetall(l2011h . 

To understand the contribution to the heating of the hot coro- 
nal plasma, interpreting the power law is crucial. Following 
the observation-based results, the power law has a negative ex- 
ponent: small flares are more numerous than large flares. If 
the frequency of small flares is high enough (the power law 
slope being steeper than -2), the heating would be dominated 
by small even ts, the nanoflares. These have been proposed by 
iParked (Il988h based on theoretical argument s and are not di- 
rectly observable with current instrumentation. IParkerldl 988h es- 
timated an energy of roughly 10'^ J (10^"* erg) and named those 
nanoflares because they are smaller than the convention of a mi- 
croflare. The latter have energies on the order of 10^*^ J and thus 
about 10"^ times the energy of a large flare. To understand how 
much nanoflares contribute to the heating or if larger flares are 
composed of smaller flares, it is important to know the distribu- 
tion of the energies of these transient events down to the lowest 
energy levels. 

The observational approach is hereby limited by a lower 
threshold given by the detection limit of the observation. In addi- 
tion, the flare energy is only a derived quantity inferred from the 
flux of exteme UV and/or X-ray photons, and is therefore only an 
indicator for the actual energy input. Energy is also lost through 
other processes, such as heat conduction or particle accelera- 
tion, and the relative merit of this energy loss mechanism might 
vary with the total energy of the respective event. Therefore even 
the power law index for the distribution of the actual energy in- 
put and for the observed energy loss might differ (for one given 
channel, e.g., the X-ray flux). 
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An alternative approach is the numerical modeling of these 
flare statistics. They are based on the theoretical assumptions 
leading to a power law distribution. A famous model uses 
the assumption that the syst em, e.g., the solar co rona, is in 
a self-organized critical state (iLu & Hamilton|[T99Tl) . The con- 
cept of the self-organized criticality was introduced bv lBak et al.l 
(Il987h . Small changes to the system will trigger further reactions 
that can be compared to avalanches (Charbonneau et al. 20()f|). 
But the question whet her the coronal magnetic field is in such a 
state remains open 

Another numerical appr oach has been employed recently by 
iDimitropoulou et all (1201 lb . They used the cellular automaton 
model to produce an average power law index of -1.80 for the 
peak energies. For that study they used several different ob- 
served active regions. However, in their model the magnetic field 
boundary remains constant over the simulation time, and the 
driving is done artificially, but both are not realistic. 

In this study we present for the first time the frequency 
distribution of transient heating events in a forward model ap- 
proach that is based on a 3D MHD model of a solar ac- 
tive region. The hot corona in this model is sustained by 
Ohmic dissipation of currents, which result from the braiding 
of magnetic field lines by convective motions prescribed at the 
lower boundary. Since the fi r st of th is type of forward model 
(iGudiksen & Nordlundll200l l2005allbh . it has been shown that 
these models not only produce a hot loop-dominated corona, 
but that they also r natch the (averaged) pro perties of spectral 
EUV observations dPeter et al.l 120041 l2006l) . Furthermore, the 
synthetic ob servations show a comparable temporal variability 
(lPeteiil2007h . and new suggestions have been made on the na- 
ture of the observ ed net Doppler shifts (?) in the transition re - 
gion and corona (iHansteen et al.l 120101 : IZacharias et al. I I20TTI) . 
These 3D MHD coronal models show that the heating is con- 
centrated in field-aligned cur rent structures and that they feature 
a strong temporal variability (iBingert & Peterl201 lb . The spatio- 
temporal distribution of the heat input also gives rise to individ- 
ual coronal loops that appear in (synthes ized) extreme UV emis - 
sion as having a constant cross section ( lPeter&Bingertll2012h . 
just as in observations. 

The success of previous models provides evidence that the 
distribution of the heating rate in space and time (on the resolved 
scales) is close to what we find on the Sun. Therefore it is timely 
to investigate details of the statistics of the energy deposition in 
the 3D MHD models, in particular the energy distribution of the 
heating events, which is the goal of this study. 

In Sect. |2] we briefly describe the 3D MHD coronal model 
and show in Sect. |3] the transient nature of the self-consistent 
coronal heat input. In Sect.|4]we derive the energy distribution of 
the heating events and discuss in Sect.|5]some observational con- 
sequences, before we investigate some aspects of self-organized 
criticality in our system in Sect|6] 

2. Coronal model 

To understand the heating of the corona, it is crucial to study 
its spatial and temporal evolution. Observations only provide in- 
formation about consequences but not directly about the heating 
mechanism and its energy deposition. Therefore we employ a 
3D MHD model of the solar atmosphere in order to analyze the 
spatial and temporal distribution of the heat input. The data used 
for the present paper is based on a coronal model introduced in 
[Bingert & Peter (201 1). The model produces a hot corona above 
a small active region with self-consistent heating. The heating in 
the model is due to Ohmic dissipation of free magnetic energy 




horizontal x [Mm] 

Fig. 1. Map of the vertical magnetic field component in the 
model photosphere at the lower boundary. Marked are the po- 
sitions of the events shown in Figs. |2] and |4] 



created by the braiding of magnetic field hnes. This mechanism 
of (DC) heating has already been proposed by iParkeii (Il983l) . 
where the braiding of field lines is a simplified picture for the 
energy transport due to the Poynting flux and the subsequent 
dissipation process. The magnetic fields anchored in the pho- 
tosphere are shuffled around by photospheric (granular) motions 
leading to a build-up of currents in the upper atmosphere. The 
plasma is then heated by dissipation of these currents. Basic pa- 
rameters derived from the numerical model, such as the energy 
flux of roughly 100 Wm"- into the corona, match observation- 
ally derived values. 

The model describes the solar atmosphere above an ob- 
served active region depicted in FigH] The domain spans over 
50x50 Mm- in the horizontal direction and 30 Mm in the ver- 
tical direction. The initial condition for the plasma parameters 
are given by a hydrostatic plane-parallel atmosphere similar to 
solar averages. The magnetic field is initially computed using a 
potential field extrapolation. 

The system is then driven by granular motions having the 
same statistical properties as the photospheric velocities ob- 
served on the sun. The evolution of the model atmosphere fol- 
lows the MHD equations. The lower boundary is prescribed 
by the granular motions and a fixed plasma pressure. The top 
boundary is hydrostatic with no heat flux. The magnetic field is 
extrapolated by matching a potential field at the top. In the hori- 
zontal directions the domain is fully periodic. 

We ran the model for one solar hour before taking any data 
to make the result independent of the initial condition. After that 
we took snapshots with a cadence of 30 seconds for a bit more 
than a hour 

For the numerical mode l we employ the Pencil code 
dBrandenburg & DobleJl2002h . It solves the MHD equations for 
a compressible, magnetized, viscous fluid. These equations are 
the conservation of mass, the momentum balance equation, the 
energy equation, and the induction equation. The Ohmic heating 
term is given as rifioj^ with j = V x BljiQ. The time-dependent 
ener gy balance in t he model includes anisotropic h eat conduc- 
tion (Spitzer '1962|), optically thin radiative losses (ICook et al.l 
1989), together with the viscous and Ohmic heating. No addi- 
tional arbitrary heating is applied in the model corona, which is 
sustained almost exclusively through the Ohmic dissipation. 
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Fig. 2. Temporal evolution of thie Oiimic lieat- 
■■B ing rate (~ j^) at four different locations (sin- 
gle grid-points) in tfie model corona: (a) grid- 
point with smallest time-averaged heat input, 
and (b) with strongest time-averaged heat input 
in the coronal volume at r>10^ K; (c) location 
at the top of a bright loop seen in the synthe- 
sized coronal emission, and (d) at the coronal 
base at the footpoint of the bright loop. The 
crosses mark the times of the snapshots. The 
geometrical heights of these four locations are 
indicated in Fig. [3] The locations are also indi- 
cated in Fig.[T]relative to the photospheric mag- 
netic field. 



3. Sample events: nanof lares and nanoflare storms 

The ohmic heating in the model is transient in time and space. 
To illustrate this, we show in Fig.|2]the temporal variation in the 
volumetric heating rate at four grid points of the numerical sim- 
ulation. These have been selected using the following criteria: 
(a) the lowest heating rate in the box, (b) the highest heating rate 
in the coronal volume, i.e., in the region where the temperature 
exceeds 10* K, (c) at the apex of a loop clearly visible in (syn- 
thesized) coronal emission, and (c) at the coronal footpoint of 
that loop. Obviously these locations span a wide range of (peak) 
heating rates, covering four orders of magnitude. 

In part, this large difference in peak heating is due to the 
strong drop in the volumetric heating rate with height. In the 
coronal part, the horizontally averaged volumetric heating rate 
drops roughly exponentially with a scale height of about 4 Mm, 
as can be seen from Fig.[3] (the heating per panicl e shows a peak 
in the transition region; cf. iBingert & Peterll20I II) . For example, 
this drop accounts for a difference in heating rate between lo- 
cations (a) and (b) in Fig.|2]of about a factor of 100 (cf. dotted 
lines in Fig. [3J. On top of this average drop with height, there 
is a strong spatial inhomogeneity caused by the structure of the 
magnetic field and how the magnetic fi eld Unes are braided. T he 
magnetic field geometry above 4 Mm (IB insert & Peterll201 1 1) is 
dominated by a large loop system connecting the main polarities 
(cf. Fig.[T]). At heights above 4 Mm, almost all field lines con- 
nect to these polarities. As a result the selected positions do not 
show any distinguished magnetic configuration. 

Most important is that the temporal variability at each posi- 
tion shows a different nature. Referring to Fig. |2] e.g., at loca- 
tion (a) the heating is due to short pulses with a lifetime of a 
few minutes, whereas at position (b) only one pulse over more 
than half an hour seems to occur At positions (c) and (d), the 
heating is composed of short pulses and of pulses with lifetimes 
of ten or more minutes. These four profiles represent the com- 
plete domain in the sense that the energy is inserted in more or 
less randomly distributed pulses. Thus it is quite challenging to 
investigate each of the million gridpoints in the numerical do- 
main. Thus one has to employ a statistical approach. Before we 
do this in Sect. H) we investigate one particular location in more 
detail. 

As an illustration we select one spatial location at the very 
base of the corona and transition region, where the (average) 
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Fig. 3. Horizontally and temporally averaged Ohmic heating rate 
(~ j~). The marked postions are at about (a) 28 Mm, (b) 10 Mm, 
(c) 22 Mm, and (d) 17 Mm according to the locations depicted in 
Fig.|2] 

heating rate per particle is highest (see iBingert & Peteijl20I ll) . 
This is shown in Fig.H] Looking at a spot where the heating rate 
per particle is at its peak gives some confidence that this is where 
most of the action is, and thus this example can be considered a 
representative one. 

To estimate the energy deposited in one single burst of heat- 
ing, the volume covered by this event is needed. The current 
structures aligned with the magnetic field are typically limited by 
the spatial resolution of the numerical experiment. Using a high- 
order scheme to solve the MHD equations, these small struc- 
tures will have an extent of roughly 4x4x4 grid points (see also 
Bingert & Peter 2011, their Fig. 9). This corresponds to a vol- 
ume of very roughly 2 Mm^ and should be considered as a lower 
limit. 

The temporal evolution in Fig. |4] shows two events in close 
succession: event (A) consists of several heating pulses occur- 
ring shortly one after another so that they merge into one sin- 
gle longer duration event lasting almost 20 minutes. The second 
event (B) consists of a single heating pulse with a larger ampU- 
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Fig. 4. Ohmic heating rate (~ j^) at a fixed point in space for more 
than a solar hour The data points are marked by crosses. Here 
a gridpoint at the very base of the corona and transition region 
(^3 Mm height) was chosen which is representative of the region 
where the heating rate per particle is maximal. The large cross 
in Fig. [2] indicates the horizontal location. The shaded regions 

(A) and (B) indicate a nanoflare storm and a single nanoflare, 
see Sect. [3] 

tude lasting only a few minutes. (The same is true for the weaker 
events at times 15 minutes and 22 minutes). 

When integrating in time and space about 10'^ J are de- 
posited during event (A) and a couple of 10'^ J are deposited 
during event (B). These energies coincide roughly with the en- 
ergy of a typical nanoflare of 10'^ J (10^"* erg) as origina lly sug- 
gested from theoretical considerations by iParkerl (Il988l ). Event 

(B) could then be "classified" as a regular nanoflare. In con- 
trast, event (A) could be related to a nanoflare storm, a phe- 
nomenon "with durations of severa l hundred to several thousand 
seconds' ' jViall & Klimchukl201 lb . Originally, these were mtro- 
duced heuristically to mod el the light curves o f warm loops with 
the help of ID models (e.g. lWarren et al.l2002h . even though this 
name was coined only later. 

Our 3D MHD model with self-consistent Ohmic heating 
shows that such nanoflare storms are a natural consequence of 
the braiding of magnetic field lines. Therefore our results sup- 
port the ad-hoc use of such bursts for the energy release in ID 
loop models. 

4. Statistical analysis of the heating rate in time 

To understand the coronal heating process, it is not sufficient to 
study single events alone. Instead we have to investigate statisti- 
cal properties of the heat input. For this we subdivide the com- 
putational domain into boxes of equal size V, and the time span 
of the simulation into time intervals of equal length r. 

As one "even?" we define the total energy input through 
Ohmic heating in one of these boxes with volume V integrated 
over time t. For example, if we choose V - (5x5x5) Mm"* and 
T = 5 min, we have 10x10x6 = 600 boxes in the 50x50x30 Mm^^ 
domain spread over 12 bins in time for a computational time of 
60 min. Consequently we will get 600x12 = 7200 individual 
events for this particular combination of V and r. We then inves- 
tigate the statistical properties of these events, in particular their 
energy distribution. We do this for a variety of combinations of 



V and T to show that our results do not depend on a particu- 
lar choice of V and r. We choose the subvolumes 0.29 Mm"*, 
0.97 Mm^, 2.31 Mm\ and 4.51 Mm\ which corresponds to 2^, 
3^, 4^, and 5^ grid points. The grid spacings are dx=d3'=0.39Mm 
and dz=0.24Mm. In this approach we use equal-sized noncubic 
volumes rather than searching for the spatial extent of each indi- 
vidual energy deposition. In Sect. |6] we provide some more de- 
tails on why this approach is useful and appropriate. A decom- 
position in real events is rather challenging and perhaps even 
impossible. The transient heating produces events that not only 
overlap in time but also overlap in space. Additionally, the en- 
ergy release shows a general decrease with height, therefore sin- 
gle events will not be easily detectable, except for a few iso- 
lated major energy releases. To reduce the investigation to those 
events, which are clearly seen in the data, would reduce the num- 
ber of events dramatically, and the statistical analysis would fail. 

4.1. Energy distribution of Oiimic dissipation in the whoie 
computational domain (from pliotospliere to corona) 

The procedure described above provides the energies of a large 
number of events. A common way to analyze this is to investi- 
gate the distribution of the energies of these events. If the fre- 
quency distribution is normalized by the event energies E, we 
get a frequency-size distribution. In addition we normalize the 
frequency-size distribution by the area and time. Here the area is 
the horizontal extent of the computational domain (50x50 Mm^), 
and the time is the duration of the simulation used in this study 
(60 min). 

Figure|5]depicts the frequency-size distributions f{E) of the 
Ohmic heating event energies for several combinations of V and 
T for the whole computational domain. The event energies range 
from 10^'^ J to 10"^ J. These are energies several orders of mag- 
nitude below and above the energy of a Parker nanoflare that 
is about 10'^ J (10^"* erg). The frequency distributions have a 
roughly constant slope of i ~ -1.2 over a wide range, irrespec- 
tive of the choice of the box size V or the temporal binning r. 
Only for high values of V and r are the number of events too 
low, and the frequency distribution is not continuous. For in- 
stance for the largest binning in space and time we have roughly 
8000 boxes, e.g. events, in comparison to the smallest binning 
with more than 10 million boxes. If the box size increases, small 
events are missing and the size of the biggest events increases. 
Therefore the covered energy range of the frequency-size distri- 
bution clearly depends on the spatial binning but depends only 
weakly on the temporal binning. 

As a consistency check, the frequency-size distribution can 
be used to derive the total energy budget of the system. The to- 
tal energy flux into the system is given by the integral of the 
frequency distribution 

f = r " f(E) E dE [WmT-] . (1) 

Assuming the frequency distribution to be a power law of the 
form /(£) = AE\ the energy flux in the range of energies from 
El to E2 is then given by 

p= _A_ £(2-)) (2) 

if 5 -2 otherwise P - A ln{E2/E[). If the slope s of the power 
law is flatter than -2, then the integral is dominated by the upper 
energy limit. Whereas for slopes steeper than -2 the result is 
dominated by the lower energy limit. 
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Table 1. Slopes of the distributions shown in Fig.|5]using maxi- 
mum likelihood estimations. See Sect. 14.21 



10-^ 

Event energy [J] 



Fig. 5. Frequency distribution of Ohmic heating event energies 
for four different spatial box volumes V and five different tempo- 
ral bin sizes t. The temporal binning t is color-coded and labeled 
in the figure. The spatial box sizes V are 0.29 Mm^, 0.97 Mm-', 
2.3 1 Mm-', and 4.5 1 Mm-'. The distribution with the smallest box 
size V is the group of lines starting at the smallest energies (to the 
left). For the two smallest binnings the number of events is large 
enough that the frequency size distribution is continuous. For the 
two larger binnings the distribution is interrupted owing to the 
low number of events. Oveiplotted is a power law with a slope 
of -1.2. The vertical dotted line ind icates the typical nanoflare 
energy as predicted bv lParkeii(ll988l ). See Sect. WJ] 



The slope of the frequency distribution of Ohmic heating 
events is -1.2 and thus flatter than -2 (c.f., Fig.|5]). When con- 
sidering the whole computational domain including the photo- 
spheric and chromospheric parts, the majority of heat is there- 
fore deposited in events with high energies. 

The energy flux into our modeled box derived from the fre- 
quency distribution in Fig.|5]based on Eq. ([T]i or (|2} is on the or- 



dO^ erg s ' cm "^). 



This is the flux at the 



derof P ^ 10*Wm-2 ( 
photospheric level, which is the lower boundary of our model 
atmosphere. This represents the averag e flux over roughly on e 
hour solar time. It matches nicely with IB insert & Peter (1201 Ih . 
where the energy flux input into the computational domain at a 
fixed time was found to be in the range between 10^ W m"~ and 
lO^Wm-2. 

It should be emphasized that these results apply to the whole 
computational domain, including the photosphere and chromo- 
sphere. In Sect. l4.3l we discuss the energy distribution and budget 
for the coronal part of the domain alone. 

4.2. Power-law test for energy distribution 

Even if the slope of the distribution functions of the Ohmic heat- 
ing events looks fairly constant over several orders of magni- 
tude, we briefly discuss how to test the validity of the assump- 
tion that this is indeed a power law. It can be problematic to 
use linear regression in a log-l og plot to retrieve the sl ope of a 
power law. Therefore we follow lVirkar & Clausetl(l20I2l) and use 
a maximum-likelihood estimation (MLE) to obtain the slope of 
a possible power law. As discussed in Virkar & Clauset (2012), 
the result is sensitive to the selection of the lower energy cut- 
off. In Table [1] we list the estimated slopes following from the 
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MLE technique for various combinations of the spatial box sizes 
V and time bins r for the events. We neglect 5% of the lowest 
energies from the analysis. 

Having derived a s lope is not yet a proof for a power- 
law. Again we follow I Virkar & Clauset (|20 1 2l) and use the 
Kolmogorov-Smimov test (lPressetal.1119921) to check whether 
the distribution functions match the power laws. The slopes re- 
trieved by the MLE give a confidence level above 0.9 for the two 
distribution functions with smallest spatial binnings (c.f, Fig.|5j. 
The highest confidence levels we find are for a slope of 1 .20. All 
slopes discussed here are around 1 .2, and the confidence levels 
are high enough to argue that the distributions of Ohmic heating 
events are consistent with power laws. 

4.3. Heating statistics for tiie coronal part alone 

We now focus on the coronal part. We define the coronal part 
either through a threshold in temperature or as the volume above 
a certain height. Both definitions (when properly chosen) give 
similar results. This will give us information about the distribu- 
tion of the heat input into the model corona. Therefore, we repeat 
the analysis from Sect. 14.11 in the subvolume of the computa- 
tional domain representing the corona and the transition region. 

To define the coronal volume, in the first case we choose 
a temperature threshold of log r/[K]=3.9 to include everything 
above the chromosphere. In the second case we define the coro- 
nal volume as being at heights above 4 Mm. This corresponds 
to the height where the horizontally averaged temperature ex- 
ceeds log T/[K] = 3.9 in our model. It is also the height where 
the scale height of the horizontally averaged volumetric heat- 
ing rate increas es and becomes cons tant for the coronal part 
(see Fig. [3] and iBingert & Peteiil20I ll) . The model atmosphere 
is very dynamic so that the isosurface for a constant temperature 
is highly corrugated. The height variation of the isothermal sur- 
face of log T/[K] = 3.9 is several Mm. This is larger than to the 
chromospheric scale height, but much smaller than the coronal 
scale height. In particular, this is much smaller than the vertical 
extent of the coronal part of our computational domain. This is 
the basic reason, why the results for these two methods will be 
essentially the same. 

When now considering the coronal volume alone (using the 
above two definitions), the event energies show a distribution 
that is composed of two power laws, with a knee at about 10'^ J, 
and a slope of -1.2 below and -2.5 above the knee (see Fig.|6]l. 
The slope at small-event energies is comparable to Fig. |5] as 
found for the whole computational domain. For energies below 
lO'** J the distribution shows the typical cut off due to limited res- 
olution and statistics. Here and in the following we only show 
the results for one particular choice of the spatial box size V 
and temporal binning t to determine the distiibution of event en- 
ergies (y=0.29Mm^^, r=5min). This is only to not clutter the 
plots, the results are independent of the particular choice of V 
and T (c.f. Sec. |6]l. The distribution in Fig. |6] represents a sub- 
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Fig. 6. Distribution of Ohmic heating event energies in the coro- 
nal volume alone. The results are shown for two definitions of 
the coronal volume: (a) logr[K]>3.9 and (b) heights z>4 [Mm]. 
Results are shown only for y=0.29Mm^^, t=5 min). Overplotted 
are lines indicating power laws with slopes of -1.2 and -2.5, 
respectively. See Sect. 14.31 



sample of the complete computational domain so that the range 
of the covered event energies is smaller than for the analysis of 
the entire domain shown in Fig.|5] 

The result that the slope of the power law is flatter than -2 
below 10'^ J and steeper than -2 above has a very pointed con- 
sequence. Below the knee most of the energy is deposited at 
higher energies, i.e., close to the knee. Above the knee most of 
the deposition is at low energies, i.e., again close to the knee. 
Therefore most of the energy will be deposited close to the knee 
at 10'^ J, which coincides with the originally proposed energy 
for nanoflares. 

A more quantitive way to locate the main contribution to the 
heat input in terms of event energies is to evaluate the energy 
deposited in different energy ranges. For this, we integrate the 
frequency distribution in Fig. |6]piecewise with a fixed interval in 
A log E [7] - 0.2 around the energy £0 

P{Eo)= f{E)EAE. (3) 

JEa-hE 

This gives the total energy flux into the corona that is deposited 
in a range of event energies centered on £()■ Figure [7] shows the 
distribution of there energy fluxes P{Eq) as a function of the 
event energy £0- The maximum flux is found at an energy of 
roughly 10'^ J (10^"^ erg) just as suspected above. The integra- 
tion is performed using an interpolation of initial distribution. 
But using the original resolution in energy does not alter the re- 
sults, and only the width of the bins would be broader The error 
for the peak flux is estimated by the half bin width of initial data, 
e.g, 0.16 or a factor of roughly 1.4 

As a sanity check, we can use the result shown in Fig. [T] 
to compute the total energy input to the model corona by 
integrating over the whole range of energies. We find this 
to be 100 Wm~^, w hich is consistent with the results given 
by iBin^ert & Peteij (l2bl 1.) and consistent with estimations 
of the coronal heating requirements based on observations 
(IWithbroe & NoveJl977h . 
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energy [J] 

Fig. 7. Energy fluxes into the corona deposited at different en- 
ergies £0 of the heating events. Derived from the distribution in 
Fig.|6] The distribution is first interpolated to allow for an inte- 
gration interval of 0.2 dex in logarithmic values. See Sect. 14.31 

This result of the distribution of energy fluxes in Fig. |7] 
shows that the upper atmosphere in our model is mainly heated 
by events with an energy of a few times lO'^ J, which roughly 
matches the typical nanoflare energy. In addition to the qualita- 
tive discussion of the position of the knee in Fig.|6] this quanti- 
tative result also provides information on the range of energies 
responsible for the energy input. From Fig. Qit is clear that en- 
ergies from 10'^ to 10'^ J (lO^^* to 10^^ erg) contribute signifi- 
cantly. 

From this we can conclude that in our 3D MHD model there 
is a preferred energy for the deposition through Ohmic heating. 
This value is the same o r der of magnitude as the nanoflare en- 
ergies derived bv iParkerl (Il988l) on theoretical grounds, despite 
all the limitations in terms of spatial resolution and magnetic 
Reynolds number in our numerical experiments. This is reassur- 
ing, because our 3D MHD mo del basically d escribes the flux- 
braiding process as proposed bv lParkerl(ll988l ). 

5. Consequences for observable quantities 

The frequency-size distributions derived in the previous sec- 
tion are based on the energy deposition due to Ohmic dissipa- 
tion. However, these events with energies close to the values of 
nanoflare energies are not directly observable on the Sun. For an 
observational study of the region where the plasma is heated, the 
only source of information we have is the emission in exteme UV 
and X-rays. And this provides only partial information on the en- 
ergy budget, because heat conduction, particle acceleration, and 
maybe other processes will also redistribute the energy, which 
will then be, e.g., radiated in other places. 

To study how the distribution of energies changes if one in- 
vestigates the brightenings related to the heating events, we syn- 
thesize the expected optically thin coronal and transition region 
emission fr om the MHD model. H ere we follow the procedure 
outlined bv lPeter et al.l(l2004l2006h . To derive reliable estimates 
of the coronal emission, the MHD model has to treat the energy 
equation appropriately; in particular it is essential to include the 
heat conduction to self-consistently set the proper coronal pres- 
sure (and thus density). Using the CHIANTI atomic data base 
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Fig. 8. Event frequency distributions of synthesized emission 
lines. The name of the hnes are given in the legend along with the 
wavelength [A] and the line formation temperature in log T/[K]. 



(??), we calculate the emission at each grid point for a num- 
ber of emission lines based on the temperature and density at 
the respective gridpoint in the model. The result is the volumet- 
ric energy loss per second through optically thin radiation in a 
given spectral line. 

We repeat the analysis of the preceding section and treat the 
emission in each of the lines in the same way as the energy input 
in Sects. l4n and 14.31 Figure [8] shows the resulting distribution 
functions for the radiative energy loss events. The energy range 
is shifted to much lower energies than in Figs. |5]and|6] This is 
mainly because each line only represents one single emission 
line, and each emission line contributes only a small fraction 
of the total energy loss. Nonetheless, the distribution functions 
show similar power laws for all emission lines with line forma- 
tion temperatures ranging from 30 000 K up to 10^ K. Owing to 
the different formation temperatures, the lines originate from dif- 
ferent densities, so that the frequency distributions have different 
cut-offs at the high-energy end. 

Most importantly, the power law shows a slope of about 
-0.8. This slope is slightly flatter at low energies and a bit 
steeper for high energies, but is flatter than the slope of -1.2 
found for the energy input. The reduction of the slope reflects 
the nonlinear connection from heat input to the final (optically 
thin) radiative output. There are basically two ways of looking 
at this phenomenon: (1) a heuristic explanation based on the en- 
ergy redistribution and (2) a more formal argument using scaling 
laws. 

At coro nal temperatu res heat conduction dominates radiative 
losses (e.g. |Priestlll982h : the higher the temperature, the more 
heat conduction dominates (below 10^ K). On average, the heat 
input (per volume and time) is smaller at higher temperatures 
at larger heights (Fig. |3]l, where the energy input is lower Thus 
heat conduction is the major loss process at the low-energy end 
of the distribution of event energies. Consequently, there is less 
radiation at the low-energy end, and the distribution of the ra- 
diative losses will be flatter, in our case the slope is only -0.8 
instead of -1.2. 

The same mechanism also explains why the knee in the en- 
ergy input (Fig. |6l Sect. 14.3b is not clearly visible in the dis- 
tribution of the radiative losses (Fig. |8). The energy deposited 



around the nanoflare energies coinciding with the location of the 
knee is mainly transported down towards the chromosphere, i.e. 
towards the high-energy end of the distribution, and is then radi- 
ated at lower temperatures in the transition region. In turn, this 
enhances the distribution at the high-energy end of the distribu- 
tion of the radiative output. 

To get a quantitative estimate of the flattening of the power 
law slope from the energy input (-1.2) to the radiative output 
(-0.8), we investig ate the scaling laws as already described by 
iRosner et al ] (ll978h . The scaling laws relate the volumetric heat 
input Eh (assumed to be constant) and the loop length L to the 
peak tempe rature T and pressu re p of a hydrostatic loop (Eq. 4.3 
and 4.4 of iRosner et"anil978h . One can easily combine these 
to relate the density n at the apex to the heat input and loop 
length, n oc e'^''l^^''. Because the dependence on L is weak, 
we discard this in the following. The energy loss s due to op- 
tically thin radiation (per volume and time) is proportional to 
the density squared, e oc and with the above scaling law we 
find s oc {Enf^''- Therefore the slope of the power law should 
change from -1.2 to -1.2x7/8 = -1.05. This change does go 
in the right direction, but does not quite explain the result of 
-0.8 we find in Fig |8l This is not very surprising, because the 
iRosner et al.l (11978 ) scaling laws apply only to hydrostatic loops 
with a constant heating rate, while our 3D MHD model is time- 
dependent with a nonconstant heating rate. 

In conclusion, the energy events for the radiative losses are 
still distributed as a power law, but the slope is flatter than the 
power law for the actual energy input. This is important for ob- 
servational studies investigating the energy input by analyzing 
the energy loss through brightenings seen in EUV and X-rays. 

6. Self-organized criticality and fractal dimension 

In Sect. |4] we derived the frequency-size distributions for Ohmic 
heating events. These distributions follow a power law over sev- 
eral orders of magnitude (c.f. Fig. |5). This power law does not 
change if the binning of the data cube is changed, either spatially 
or temporally. This behavior is typical of a power-law frequency- 
size distribution, so it indicates that our modeled energy input 
also follows a power law. Systems that are characte rized by a 
powe r law are called self-organized critical systems (iBak et al.l 
[l987l) . Self-organized criticality can also be applied to other sys- 
tems such as 1 // noise, and it shows the scale invariance of the 
system on a macroscopic level. This scale invariance manifests 
itself by the fact that the slopes of the power laws are indepen- 
dent of the spatial scale range under consideration, i.e. of the 
spatial binning. Therefore the approach of using a constant bin- 
ning for the analysis in Sects. 14.11 and 14.31 does not affect the 
results. The constant binning refers to equally sized volumes. 
Likewise the noncubic boxes do not alter the result since we can 
use arbitrary volumes in a scale-independent system. 

We find in our analysis that the frequency size distribu- 
tion of the event energies can be represented by a single slope 
alone. This is cons i stent w ith other model approaches , e.g, 
iDimitropoulou et all ( 201l]) and Georgoulis et al. (1200 ih . and 



with observations ("Aschwanden & Freeland' '2012'). In earlier 
studies lOeorgoulis & Vlahos (1996) and Georgoulis & Vlahoi 
(Il998h used cellular automata models and find a double power 
law behavior for the event size. That our distribution functions 
follow a single power law can be accounted for by the short time 
period of one hour and by the small field of few of our model 
active region. 

Following lAschwandenl ( 1201 Ol) the slope of the power law 
of the distribution functions point to the propagation properties 
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of the instability leading to the dissipative events. The fractal 
dime nsion D is thereby correlated to the slope 5 by D - 3/s 
dAsc hwanden 2010; Aschwanden & Freeland 20121). Since we 
measure a slope of 1 .2, the fractal dimension of the system would 
be 2.5 for the Ohmic dissipation events. This means that the frac- 
tal dimension is not bound to the 2D surface of separatrix layers, 
and that the propagation of the instability cannot spread freely in 
all three dimensions. 

Formally, using the slope of -1 of the power law for the 
radiative losses (c.f. Fig. |8]i, we find a fractal dimension of 3. 
However, because the emission is just a secondary quantity fol- 
lowing from the heat input and strongly affected by the heat con- 
duction, the diagnostic value of this result is limited, though. 



7. Conclusions 

In Sec. Owe discussed the transient nature of the heating based 
on time series of the Ohmic heating at specific locations. There 
the temporal evolution of the heating rate shows a stochastic 
nature, as well as a qualitative difference at different places. 
We saw long bursts (10 min and longer), which may consist 
of several short pulses. At other locations we saw only short 
pulses (1 min and less). These events could correspond to the 
nanoflare storms and nanoflares often assumed in ID loop mod- 
els ( Viall & Klimchuk 201 1). We used these time series only for 
a qualitative description and for a rough estimate of the magni- 
tude of heating events. We saw that these events deposit energies 
on the the or der of 10^^ J, w hich agrees with the nanoflare energy 
estimated bvlParked (119881) . 

We also presented a more sophisticated approach using sta- 
tistical methods to analyze these events in the complete domain. 
The statistical analysis of the coronal heating reveals a power 
law for the heating event energies. Even though the driver, i.e. 
the granular motion, does not work on all scales, the slope of the 
power law is constant over several orders of magnitude and flat- 
ter than -2. We find that in the transition region and corona the 
dominant energy for the heating events is the nanoflare energy, 
i.e. 10'^ J. Since the higher energy events occur on average at 
lower heights, the heating is footpoint-dominated. The event en- 
ergies we find in our model are below the detection limit of cur- 
rent instrumentatio n and support the scenario of nanoflare heat- 
ing as proposed bv lParkerl(ll988b . 

In our model the energy provided by the field-line braiding 
is converted into heating through Ohmic dissipation and not into 
particle acceleration. This implies that non local effects are not 
incorporated. This is justified by the electron mean-free-path be- 
ing smaller than the grid spacing of the model. The energy is 
redistributed by advection and conduction or radiated by opti- 
cally thin emission. From our model atmosphere we synthesized 
observable extreme UV emission lines. The statistical analysis 
of these lines shows again a power-law behavior The slope of 
this power law is around 1 and smaller than the slope of the 
heating process. Even though the power law is maintained over 
several orders of magnitude, in contrast to the heating process, 
the radiative losses are not in a self-organized critical state. This 
is because the radiative losses are coupled to the heating in a 
nonlinearly. 

In this study we showed that a 3D MHD model of the ac- 
tive region solar corona has a preferred energy for depositing 
of heat through Ohmic dissipation following field-line braiding. 
The preferred energy of about 10 J (1 0^"^ erg) corresponds well 
to the value derived by iParkerl (1 19881) and thus is yet another 
piece of evidence that, despite all their shortcomings, the 3D 



MHD models for the flux-braiding mechanism provide a good 
description of the heating of the solar corona in active regions. 
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